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Abstract 

A dressed state approach to mixing of bosonic matter waves is presented. Two cases are stud- 
ied using this formalism. In the first, two macroscopically populated modes of atoms (two-wave 
mixing) are coupled through the presence of light. In the second case, three modes of Bogoliubov 
quasiparticles (three-wave mixing) are coupled through s-wave interaction. In both cases wave 
mixing induces oscillations in the population of the different modes that decay due to interactions. 
Analytic expressions for the dressed basis spectrum and the evolution of the mode populations in 
time are derived both for resonant mixing and non-resonant mixing. Oscillations in the population 
of a given mode are shown to lead to a splitting in the decay spectrum of that mode, in analogy to 
the optical Autler-Townes splitting in the decay spectrum of a strongly driven atom. These effects 
cannot be described by a mean-field approximation. 

PACS numbers: 03.75.Gg 
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I. INTRODUCTION 



Wave mixing is a well known phenomenon which occurs in nonlinear systems. The realiza- 
tion of Bose-Einstein condensates (BEG) in atomic vapor paved the path to the experimental 
study of mixing of matter waves [l,Q, [J, in which the s-wave interactions are the cause for 
the required nonlinearity. This exciting new effect is analogous to the wave mixing of opti- 
cal modes in nonlinear processes such as parametric down- conversion |4|. As in the optical 
case, atomic wave mixing is predicted to generate number squeezed states. These states 
are usually referred to as quantum states with a well defined number and a spread in the 
phase, contrary to the BEC itself, which is described well by a classical wave in the sense 
that both number and phase are well defined. These squeezed states can be used to perform 
sub-shotnoise measurements in interferometry experiments jsj. 

In previous work, the nonlinear mixing of Bogoliubov quasiparticles was studied using 
the Gross-Pitaevskii equation (GPE) 6]. In the high momentum limit, in which most 
experiments are performed, atomic four- wave mixing was studied using the GPE [7|, and 
analyzed directly in the quantum many body formalism using angular momentum for mixing 
of different hyperfine modes . The coherence in wave mixing of two light modes with two 
atomic modes was studied experimentally £ | and demonstrated when only one atomic mode 
and one light mode are initially populated lol. \v\. 

Many dynamical effects in ultra-cold atoms can be described by the time dependent GPE, 
which governs the evolution of a macroscopically occupied wavefunction. This approach 
cannot be used in order to study quantum effects where number and phase cannot both be 
determined at once. Here we develop a dressed state model [l^j], which unlike the GPE, does 
not assume all atoms are in the same single particle state. We use this model to describe 
effects that are beyond mean-field, such as dephasing due to quantum uncertainty in Rabi 
oscillations between momentum modes, and the decay of a matter wave due to interaction 
with the quasi-continuum of unoccupied modes. 

Weak interactions are usually taken into account by the Bogoliubov transformation from 
the atomic basis to a basis of quasiparticles. Interactions between quasiparticles and their 
decay were first analyzed by Beliaev jl3j |. and experimentally verified 14, 15]. In this paper 
we study wave mixing between macroscopically occupied quasiparticle modes. For simplicity 
we consider throughout the paper only homogeneous condensates at T = 0. We find that 
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interactions between quasiparticles can lead to new phenomena in matter waves, analogous 
to phenomena studied in connection to interactions of light waves with matter, such as a 
splitting in the decay spectrum, power broadening and energy shifts [16]. We first apply the 
dressed state formalism to the simple case of two matter waves coupled by laser fields. We 
then adapt it to solve the more interesting case of three-wave mixing as a locally perturbed 
two-wave mixing. 

Treating the wave mixing problem in the Bogoliubov basis allows us to consider mixing 
between quasiparticles in the low momentum regime, where the basis of atomic modes (used 
e.g. to describe four- wave mixing of matter waves [l|, |3j) is not appropriate due to the large 
— k atomic component of a quasiparticle with momentum k. Such a low momentum regime 
is of particular interest since both the density of states and quantum interference effects yield 
very low decay rates ^3, 17], and is therefore more suitable experimentally for the study of 



relative number squeezing and entanglement between matter waves. In the high momentum 
regime, three-wave mixing coincides with the conventional four-wave mixing of matter waves 
as long as the condensate, which is the fourth wave, is not substantially depleted. 

The layout of this paper is as follows: First, in section |H] we review the Bogoliubov 
transformation and the wave mixing Hamiltonians for two physical scenarios. We describe 
both the mixing of Bogoliubov quasiparticles due to the s-wave interaction, and the mixing 
of atomic modes which are coupled by laser fields, and in which case interactions may be 
included as a perturbation. We also derive the perturbative dynamics and energy shifts due 
to the interactions between Bogoliubov quasiparticles. Next, in section we solve the 
two-wave mixing Hamiltonian. For noninteracting atoms, this is in fact a single particle 
problem. However, we use the second-quantized formalism, which will become a necessity 
once interactions are introduced, and nonlinear dynamics emerge. We compare this formal- 
ism to the two- mode Gross Pitaevskii model Q], and discuss the dynamical instability, as 
well as a decay due to quantum uncertainty in the number of quasiparticles. In section 
IIVI we turn to the nonlinear three-wave mixing between Bogoliubov quasiparticles, using 
the framework introduced in section IIII1 We study the dynamics and spectrum of resonant 
and non-resonant three-wave mixing. We find an analytic approximation which is in good 
agreement with numerical diagonalization, and derive the scaling laws which let us describe 
a realistic system by matrices which can be diagonalized on a personal computer. Finally, 
in section El we discuss the decay of a quasiparticle coupled both to a quasi- continuum and 
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an additional macroscopically occupied quasiparticle mode. We find a splitting in the decay 
spectrum in analogy to the Autler-Townes splitting of the spectrum of an atom undergoing 
Rabi oscillations. We conclude with a brief summary and examples of other cases to which 
our formalism can be applied. 



II. WAVE MIXING HAMILTONIANS FOR ULTRA-COLD ATOMS 

In this paper we concentrate on a homogeneous system of ultra-cold bosons, with s-wave 
interactions. We use this system to derive the physical Hamiltonians giving rise to two and 
three-wave mixing. The three- wave mixing is of Bogoliubov quasiparticles over a condensate 
which are inherently coupled by interactions. It therefore describes interactions correctly 
even for small momentum modes, but is not applicable to strong excitations which deplete 
the condensate 19]. The two- wave mixing is of two macroscopically populated momentum 
modes, which can be depleted, and are coupled by an external laser field. Interactions 
are incorporated as a mean- field shift, thus the two- wave model is not applicable for small 
momenta. 



A. Wave mixing of Bogoliubov quasiparticles 

In momentum representation the Hamiltonian governing a system of ultra-cold interacting 
bosons is given by 20]: 

H = J2 + 2V ^ a k«l«mak+l-m. (1) 

k k,l,m 

g = A7fh 2 a s /M is a constant proportional to the scattering length a s , describing the atom- 
atom interaction. V is the volume of the condensate and M is the atomic mass. Assuming 
a macroscopic occupation of the ground state N , thus replacing both a and at by V^o 



this Hamiltonian is diagonalized to quadratic order by the Bogoliubov transformation 

flk = u k b k - v k P_ k (2) 
a{ = u k b\. - Wfcfo-k- (3) 
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When taken to cubic terms in b, Eq. ((TJ can be written as H = Hq + Hi nt with 

H = E g + J2tk>blb h ,, (4) 

k' 

where E g is the energy of the Bogoliubov vacuum, and 

Him = ^77^ [ Ak '^' (^'^q'^k'-q' + &£'&q'&k'-q') + #k',q' (^'^^-(k'+q') + &k'M-(k'+q')) 

k',q' 

(5) 

For weakly excited condensates H int is typically small, and is usually neglected, as in the 
derivation of the Bogoliubov spectrum = \J ~^f(J§j^ + 2/i), where // = gp is the chemical 
potential and p is the condensate density. The quasiparticle modes which diagonalize H 
interact via H int . The terms ^4k,q ; -Bk,q are combinations of the it's and v's originating from 
interference between different atomic collision pathways, and are given by 

-4k,q = 2u k (u q U\^-q\ - MgW| k -q| ~ fgW|k-q[) ~ 2l>fc (fgf |k-q| - V^k-ql ~ « ? f|k-q|) (6) 
-B k ,q = 2u k (f ? f|k-q| - « 9 f[k-q| ~ UgW[k-q[) ~ 2^fc (« 9 W|k-q| ~ Ug«[k-q| ~ ^|k-q|) • (7) 

ylk, q and i?k,q are symmetric under the exchange q «-> k — q, as implied by Eq. (jHJ) . The 
different terms of H int describe different processes resulting from the quasiparticle interac- 
tions. These processes are described diagrammatically in Fig. ^ The first term in Eq. (J5J) 
corresponds to the decay of one quasiparticle into two, and is analogous to parametric down 
conversion. This decay is referred to as Beliaev damping. The second term describes two 
quasiparticles colliding, and forming one quasiparticle with the total momentum. This term 
governs the so called Landau damping. The last two terms in Hi nt (proportional to -Bk',q') 
describe the spontaneous formation and annihilation of three quasiparticles with zero total 
momentum. 



1. Resonant wave mixing 

Consider the case where there are initially iV « iVj excitations in the momentum 
mode k, and no other excitations. The state |iVk) is degenerate under Hq with the states 
\(N — l)k, l q ', lk-q') for those q' that fulfill the resonance condition 

e ? ' + e|k-q'| = efc- (8) 
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FIG. 1: The four processes described by Hint in Eq. ©: (a) Beliaev damping (b) Landau damping 
(c) spontaneous creation of 3 quasiparticles (d) spontaneous annihilation of 3 quasiparticles. Only 
processes (a) and (b) can be resonant. 

The mode k can therefore decay via the Beliaev term in Hi nt into pairs of modes on a shell 



The factors of 2 in Eq. Q are due to bosonic exchange symmetry between q' and k — q'. 
The FGR is valid as long as the quasi-continuum contains enough modes N c such that 
during the time of the experiment the average population of these modes will be less than 
1, i.e. Tt <C iV c . The second (Landau) term in Hi nt does not contribute to the evolution 
at temperatures much lower than the chemical potential, since there is only a negligible 
occupation of modes other than k. The third and forth terms have no contribution to the 
decay since they do not conserve energy. 

Now consider the evolution of the same state |iVk), when there is an initial seed of M 
excitations of mode q which satisfies the resonance condition [Eq.flSJ], i.e. it is on the 
energy-momentum conserving shell. Our initial condition is thus |(A r )k, (M) q ). Such a 
situation is presented schematically in momentum space in Fig. EH- The condensate, which 
is the quasiparticle vacuum is represented by a dark gray ellipse. The initially occupied 
quasiparticle modes are plotted as light gray ellipses, and the initially unpopulated mode 
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k — q is represented as an empty ellipse. The line describes the quasi-continuum of possible 
modes for Beliaev damping from mode k. For M N c the initially occupied mode q will 
merely enhance the decay rate in Eq. © to 

„ 2nMg 2 N NA kQ , . 

T " = T + ^r-^w^- (10) 

This enhancement of F is analogous to the optical power-broadening, where an occupied 
mode of the electromagnetic field broadens the natural linewidth of an atom T, which is also 
the decay rate due to spontaneous eruissiou Q. If, however M > S n a situatiou that is 
experimentally achievable in cold atoms where N c ~ 10 3 — 10 4 , the dynamics are changed 
significantly from the simple decay. The term containing k, q in the sum will dominate H int 
in Eq. (J3J) and decay into two unpopulated modes can be neglected, yielding H int ps H^, 
where we define the resonant three-wave Hamiltonian as 



Tires 

*3W- 2 y 



^(^A+^i^)]. (ii) 



Since the modes q and k — q will both become macroscopically occupied, the second term 
in Eq. (fTT]) cannot be ignored. This case will be studied in section IIV1 Since the states 
| (N — n)k, (M + n) q , nk-q) with different n are degenerate under H , the eigenstates of 
in this subspace are also eigenstates of H, and once calculated, both dynamics and spectra 
are extracted in a straightforward manner. In sections ll VI and IVl H^ v will lead to nonlinear 
oscillations, and a splitting in the decay spectrum of the mode k. 



2. Non-resonant wave mixing 



The non-resonant terms of Hi nt in Eq. (JHJ) do not contribute to the perturbative dynamics. 
They do, however, lead to energy shifts. Using second order perturbation theory we find an 
energy shift of the Bogoliubov quasiparticle energy of mode k by 
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(12) 



where p stands for the principal part. This energy shift is analogous to the Lamb shift of 
atomic states due to the Electromagnetic vacuum ^(|. In analogy to the resonant wave 
mixing, we study the Hamiltonian (|5|). again with the initial condition |(iV)k, (M) q ). This 
time we consider the case where the mode q is not on the energy-momentum conserving 
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FIG. 2: A schematic momentum-space description of three-wave mixing, (a) The M-fold occupa- 
tion of the seed mode q fulfills the resonance condition of Eq. (jSJ , and is therefore on the solid line 
which marks the quasi-continuum of possible states for Beliaev damping of mode k. (b) The seed 
mode q is not on the energy conserving shell, leading to non-resonant wave mixing. In both figures 
the condensate is marked as a dark gray ellipse, the initially occupied modes are marked as light 
gray ellipses, and the empty modes which fulfill momentum conservation are represented as empty 
ellipses. The shell of possible modes is calculated for k£ = 3.2, and is therefore not completely 
spherical as expected for atomic s-wave scattering, but rather lemon-shaped Note that only 
in the high momentum limit do the plotted modes correspond to atomic momentum modes. 

shell defined by Eq. (jSJ), but is rather detuned by HA = e q + £|k-q| _ ^k, as schematically 
shown in Fig. We consider small values of detuning HA << e&, e q , hence we neglect the 
terms -Bk,q m Eq. (jSj) since the energy denominators of the -Bk,q terms are much larger than 
those of the Ak,q terms. This is in analogy to the rotating wave approximation in optics. 
To keep the states \(N — n) k , n k ^ q , (M + n) q ) degenerate under H , Eq. (jlj) is modified to 
the form H Q = E g + Y^w e fc'^k'^k' — ^jk- q |^l k -q|^- ^ n or der for H to remain the same as in 
Eq. (jSJ) we must generalize H^w to 

Hw = ^ [A k ,q + SkSt gt_ q )] + Sf k _ q| 6 |k _ q| A. (13) 

We will show in section IS that Hi nt leads to a splitting in the damping spectrum of mode 
k. The detuning leads to an asymmetry in the peaks. When the detuning is large (but still 
smaller than e^/fr) , one peak vanishes, and the other peak dominates. The large peak is 
shifted from due to the macroscopically occupied mode q by 

g*N Q M |A k , q | 2 

6€k = 2A- (M) 

This shift in the spectrum is analogous to the ac-Stark shift and appears in addition to the 
Lamb shift of Eq. Q. 
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B. Wave mixing due to Bragg coupling 



The interaction Hamiltonian described in Eq. (jSJ) takes into account the quantum de- 
pletion of atoms from the k = mode in the many-body ground state of the BEC due to 
interactions. In fact, this effect is what leads to the momentum dependence of A^. How- 
ever Eq. (J3J) is derived under the assumption that at all times the population in all modes 
much is smaller than Nq, i.e. the condensate is undepleted. There are cases where the 
condensate is depleted, for example by a strong two-photon Bragg coupling leading to Rabi 
oscillations between momentum modes and ke [19]. The Bragg coupling is characterized 
by a frequency ujb and a wavevector kg, that are the difference in the laser frequencies and 
wavevectors respectively and by an effective Rabi frequency Q which describes the coupling 
strength. The detuning of the Bragg pulse is 5 = ujb — hk B /2M. In the rotating wave ap- 
proximation (5 <C Hk B /2M) the resonant Bragg pulse couples the macroscopically occupied 
mode k = to kg, so we may estimate the effect of the Bragg pulse as j^lj: 



where we have neglected coupling to all other momentum modes (—kg, 2kg, . . .) which are 
far off resonance. So far interactions are not included in this subsection, and therefore we 
use the atomic basis rather than a quasiparticle basis. As in Eq. (jlHj) . for a detuning from 
resonance that is much smaller than the detuning from other momentum modes (\S\ <C 
\wb ~ h{mkB) 2 /2M\ for m ^ 1), the Bragg Hamiltonian is simply modified to 



Here = a^a^ is the number operator of mode k. Neglecting interactions between the 
atoms, H 2 w describes the single particle two level system. If one adds H 2 w to the Hamil- 
tonian (JTJ), it is no longer possible to use the Bogoliubov approximation, since the Bragg 
process depletes the condensate. One can however include atomic interactions by taking 
only the k = and the k = ke modes into account in the interaction term in Eq. This 
amounts to neglecting terms proportional to v\ in the Hamiltonian, such as the quantum de- 
pletion of the condensate. Neglecting Vk and approximating ~ 1 is a good approximation 
for momenta larger than the inverse of the healing length £ = (8n pa s )~ 1 ' ' 2 . In this regime the 
Bogoliubov quasiparticle energy is a free particle parabola shifted by the mean-field energy 




(15) 




(16) 
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fi, and the quasiparticle wavefunction is equal to the free particle wavefunction |20|. The 
coupling Hamiltonian including the resonant Bragg coupling and the atomic interactions of 
momentum modes and ke can be written as: 

H 2W = ~Y («o«k s + ao4 B ) + 7^7 ["o («o - 1) + ™k B ( n kg - 1) + 4n n ks ] • (17) 

Using the conservation of particles no + ^k s — N, the interaction term in Eq. (|17|) can be 
simplified to (A^ 2 — N + 2fiQn^ B )g /2V . For a small excitation, i.e. n-^ B <C no, this coincides 
with the mean-field energy nk s /i of Bogoliubov quasiparticles, while for a nearly depleted 
condensate, where n <C n^ B1 the energy is the same due to the symmetry no «-> nk B . 
This situation describes a condensate that is moving in the lab frame of reference with 
momentum k# and no excitations with momentum — kg relative to the condensate. Here 
the mean field energy per quasiparticle is gn^/V . H % ^w is no longer a linear Hamiltonian 
and cannot be solved in first quantization as a single particle problem. In fact, even the 
Gross-Pitaevskii equation is inadequate for solving this Hamiltonian for certain parameters 
Q| . The Hamiltonians (jl5H17j) are studied in section IIHI Systems governed by Eq. (JT7jl 
also exhibit oscillations, a splitting in the excitation spectrum and a variation in the decay 



spectrum. These phenomena have been demonstrated experimentally 
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III. TWO- WAVE MIXING 

In this section we discuss the diagonalization of the two-wave mixing Hamiltonians (|15l - 
H2j). We start with Eq. (fTHjl . Although this is a single atom problem (since it includes no 
interactions) we diagonalize the Hamiltonian in second quantization, in order to develop a 
notation used later for problems such as three-wave mixing which are many-body in nature. 
In subsection IIII CI we use this simple notation to describe the perturbative dynamics of 
a many-body problem including s-wave interactions between the atoms oscillating between 
large momentum modes. 

Although H 2 w involves two fields, the conserved quantity A" = aja + a^a^ allows us 
to diagonalize in a subspace with one quantum number n = nk B ,varying between and N, 
and representing a Fock state \n) = \(N — n)o, ^k B )- 

H 2 w couples between states in this subspace which have a difference of 1 in the quantum 
number n, and can be represented in the resonant case by the (A^+ 1) x (N + 1) tri-diagonal 
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FIG. 3: (a) Absolute value squared of the transformation matrix between the Fock basis and the 
dressed basis which diagonalizes ifiifv given in Eq. (|15|) |(n|ma;)| , for N = 100. Dashed line is the 
circle (n — j) 2 + m 2 = j 2 . The thin horizontal lines show the location of the n = 7 and n = N — 7 
rows which are plotted in Fig. 0J (b) The spectrum of which is strictly linear. Energy is in 

units of Ml 

matrix 
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rrres _ m 
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V 



y/(N-0)l 



(18) 



When is diagonalized, we get a new set of + 1 eigenstates, that are dressed by the 
interaction, labelled by the quantum numbers m x , which can take + 1 different values. 

The squared absolute value of the matrix elements of the transformation matrix between 
the Fock states |n) and the dressed states m x , are shown in Fig. El for A^ = 100, along with 
the corresponding eigen-energies. As can be seen, the spectrum is linear, E mx = hQNm x , 
indicating this is indeed a single particle problem. Once the transformation matrix is solved, 
dynamics can be calculated trivially by projecting the initial state on the dressed basis 
and evolving each dressed state according to \m x {t)) = e~ lEmxt ^ h \m x (0)), yielding the well 
known Rabi oscillations, with frequency Q. These oscillations have been demonstrated 
experimentally 
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A. The Schwinger boson mapping 



We now turn to describing the problem in terms of angular momentum. For non 



24] 



interacting resonant coupling the exact solution was found by Tavis and Cummings 
This approach will prove to be useful later when we deal with problems which are no longer 
linear. The two dimensional harmonic oscillator can be solved in different bases. Each oscil- 
lator can be solved separately and the eigenfunctions are tensor products of one dimensional 
oscillators. Another option is to work in polar coordinates, and find eigenstates of angular 



momentum and of the radial equation. The Schwinger boson mapping [25|, |26( maps the 
Fock basis of decoupled one dimensional harmonic oscillators, with a fixed sum of excita- 
tions N = n + rikg onto the angular momentum basis with a fixed total angular momentum 
j = N/2. The mapping is defined by the following operators 



J+ = a ai B , (19) 
it 



4 

J~ = a kj3 al, (20) 
= \ ( a L ak B - a o«o) • (21) 

Both representations fulfill the angular momentum commutation relations, and this mapping 
is completed by associating the quantum numbers 

j ee nk * +n ° = * (23) 
J 2 2 y 1 

Thus a state \(N — n)o, n^) in the Fock basis is equivalent to the state 
\j = N/2,m z = nk B — N/2) in the angular momentum basis. Using this mapping we write 
the Hamiltonian in the angular momentum basis simply as 

iz™ = + j~) = hnj x . (24) 

The eigenstates of H^, which we refer to as dressed states, are labelled by their projection 
on the x axis m x = —j, . . . , j, rather than their projection on the z axis m z = n — j, and 
the spectrum is linear E mx = Mlm x . The transformation matrix from the Fock basis to the 
dressed basis, shown in Fig. EJ is also straightforward to understand in this representation, 
and is simply a N+l dimensional representation of an SO (3) rotation around the y axis by an 
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angle 9 = — 7r/2 with matrix elements, known as the Wigner matrix elements, dP m mz (— 7r/2), 
which are defined analytically in j^. 

In the limit of large j, we can use classical reasoning to predict the shape of the trans- 
formation matrix. A classical spin of length j with a projection m z on the z axis could 



be found anywhere on a circle of radius yj 2 — m 2 . in the x-y plane with equal probability. 
Therefore, the probability to measure a certain projection m x on the x axis vanishes for 
\ m x\ > a/j 2 — rn 2 and is proportional to (j 2 — m 2 ) -1 / 2 for \m x \ < \J j 2 — m 2 . For this rea- 
son we expect that for large enough representations, the transformation matrix should be 
bounded by the circle m 2 + m 2 = j 2 , and its amplitude should peak on this circle. In Fig. 
Of a), this circle is plotted as a dashed line. Indeed, the largest amplitudes are found on this 
circle, and the amplitudes outside are completely negligible. 

i?2VK 1S analogous to the Hamiltonian of a magnetic field in the X direction applied to a 
spin j particle. The spectrum presented in Fig. Efb) is therefore analogous to the Zeeman 
energy of that particle. 

If one of the momentum modes, say 0, has a large occupation, the operator aJak B + a k ao 
can be estimated by \/no(a kB + dk B ) which is proportional to the position operator for the 
mode ke, whose eigenfunctions are the single harmonic oscillator wavef unctions. Therefore 
for n close enough to the representation of the Fock states in the dressed basis are indeed 



very similar to the eigenfunctions of a single harmonic oscillator and are given by 



28,Q 



(n\m x ) = C,n-i(V2) « (-l)"r V V (V7 arcsin {m x /j)) , (25) 

where u n (x) is the nth eigenfunction of the harmonic oscillator. Figure 0] presents the 
overlap of the Fock state n = 7 with the dressed basis m x , (n = 7\m x ). The overlap is seen 
to fit Eq. (J25J) remarkably well. Using the symmetry d J mm , = {—l) m ~ m ' d 3 _ m _ m , we find 
that (— i^ n ~ J ~ mx (jV — n\ — m x ) = (n\m x ) for all m x . This reflects the symmetry <-> ke in 

rjres 
n 2W- 



B. Detuning 

Next, we describe the case in which there is a detuning between the Bragg frequency and 
the difference in the energy of the momentum modes. We therefore add the detuning to 
H2W as we did in Eq. (|16j). In the matrix (|18|h a term —nh5 is added to the nth element 
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FIG. 4: The overlap between the Fock state n = 7 and the dressed state of H'^ m x , (n = 7\m x ) = 
(— l) n_ J _ma; (n = N — 7\ — m x ) for N = 100. The harmonic oscillator approximation given by Eq. 
(|25|) (line) is seen to be a good approximation. 
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FIG. 5: (a) Absolute value squared of the transformation matrix between the Fock basis and the 
dressed basis of H2W, \( , n\' m x)\ 2 , for N = 100 and 8 = 0,. Dashed line is the shape (n — j — 
(5/fl)m x ) 2 + m 2 = j 2 , which is the generalization of the shape in Fig. IHIb) The spectrum of H2W1 
which is again strictly linear, with slope V^ 2 + 5 2 , centered at —jh5. Energy is in units of HQ. 

on the main diagonal of the Hamiltonian. The matrix is diagonalized by the transformation 
matrix shown in Fig. along with the spectrum of the Hamiltonian ((TSj) for N = 100 and 
5 = Q. Again, the linear spectrum reflects the single particle nature of the Hamiltonian. 

In order to analytically solve Hamiltonian ()16j) we use the fact that = N/2 + J z = 
j + J z , to write H 2 w 111 the angular momentum representation as 

H 2W = MIJ X - hS (j z + j] = hVW + 5 2 Ji - j5. (26) 

The energy spectrum is again linear, with a spacing of H\/Q 2 + S 2 , centered about —j8. The 
eigenstates are those of a spin component operator pointing in the I direction in the x — z 
plane, the angle from the z axis being 9 = 7r + arctan (Q/5). Thus the transformation matrix 
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is once again a rotation by 6 about the y axis. The shape of the transformation matrix is 
again captured by classical reasoning (dashed line in FigE^). 

Dynamics are also easily understood in the angular momentum picture, in terms of spin 
precession. The Fock basis is characterized by a quantization axis pointing up on the Bloch 
sphere, while the dressed basis has the quantization axis tilted by an angle 9 towards the 
x axis. This gives a visual understanding of the transformation matrix, as well as the 
dynamics. The pseudospin evolves by precession about the axis I, with frequency y/fl 2 + S 2 . 
The well known properties of Rabi oscillations are obtained in a straight-forward manner. 
Inversion can only be achieved when the driving field is resonant to the transition between 
the two modes. As the detuning increases, the amplitude of oscillation decreases and the 
frequency, which is simply the energy difference between dressed states increases. Since we 
have neglected interactions, the spectrum is strictly linear, and there is no dephasing. 

As in the resonant case discussed in section IIII Al the symmetry of the rotation matrix 
allows us to relate the row n to N — n. For small n, both are approximated by properly 
translated and stretched harmonic oscillator wavefunctions 1281. 



C. Rabi oscillations of interacting atoms 

Upon inclusion of interactions between the atoms, the problem ceases to be a single body 
problem, and the second quantization method reviewed in the previous sections becomes 
useful. The inclusion of interactions in the resonant Bragg coupling of Eq. (|17|) is correct in 
the limit kB<£ > 1 as discussed in section III Bl In this subsection we assume for simplicity 
5 = 0. Generalization to non-resonant Bragg coupling is trivial. 

In the angular momentum representation, the interaction term in Eq. (|T7j) can be rewrit- 
ten up to a constant which we neglect in the form 

H^ = hnj x + (fi/N)J 2 z . (27) 

The transformation matrix obtained by numerically diagonalizing H^w is shown along with 
the obtained spectrum in FigElfor /i = Ml/ 2. 

The interactions are seen to change the circular shape of the transformation matrix to 
an oval shape, breaking the symmetry m x «-> —m x . The spectrum which is no longer linear 
indicates this is not a single particle problem. 
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FIG. 6: (a) Absolute value squared of the transformation matrix elements, for H^w with fj, = Ml/2 
and N = 100. The dashed line is the perturbative classical solution m x = ±y/ j 2 — m 2 + /j,m 2 /4:j. 
(b) The spectrum (solid line) is no longer linear. The non-interacting spectrum is shown for 
reference as a dashed line. The first order perturbation of Eq. (|29j) is indistinguishable from the 
solid line. Energy is in units of Ml. 

In order to describe analytically the effect interactions have on the dynamics, we use first 
order perturbation theory to calculate the energy correction to the state \m x ). We start 
with the eigen-energies of the non- interacting case Eml = H£lm x , and add the correction 

E^l = (m x \ H int K> = - (m x \ ^J 2 Z K> = £ (m x \ i - J 2 (J + 1} \m x ) (28) 

The overall energy of the state \m x ) is therefore 

E mx = H J -^- + ?ittm x + ^m 2 x . (29) 

This perturbative result for the spectrum is indistinguishable from the solid line in Fig. 
which is obtained by exact diagonalization. The dashed line in Fig. EJa) is a perturbative 
correction to that of Fig. Efa). It is derived by assuming the energy is nearly linear, 
and therefore (fi/hQ)ml/4:j is added to m x . This line is seen to indeed describe the peak 
amplitudes of the transformation matrix very well. 

The departure of the energy spectrum from a linear one, as shown in FigEb; leads to a 
dephasing between the different dressed states. According to Eq. (|'25|). an initial Fock state 
n = is a gaussian with width y/j. Hence a collapse should occur at a time scale over which a 
7r phase is acquired in the quadratic terms of the energies of the states m x = and m x = y/j. 
According to Eq. (|2*U|) this energy is /i/4, leading to a collapse time of t co u apse = 4nh/fi. In 
Fig [7[ we see the decaying oscillations due to the perturbative dynamics governed by Eq. 
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FIG. 7: The expectation value for population of mode ke evolving according to H2W mt with 
interactions for [i = HQ/2, and N = 100. The dots are according to the numerical diagonalization 
of Ml/2J X + n/NJ^. The line is according to H<iw with the first order correction to the energies 

(HI. 



fit the exact numerical solution remarkably well. A revival occurs when each quadratic 
phase is a multiple of 2tt. This happens for multiples of the t rev i va i = 4TrhN/fi. In realistic 
experimental conditions, N is typically on the order of 10 5 and thus decoherence is likely to 
occur much before t reviva i and we do not expect revivals to be seen. 

The broadening in the spectrum which causes t coUapse is proportional to the chemical 
potential, and exists for a spatially homogeneous system. It should not be confused with 
the broadening in the resonance of Bogoliubov quasiparticles over a spatial inhomogeneous 



BEC, known as inhomogeneous broadening 30J, where the local density gives rise to a local 



chemical potential and a local mean-field shift. For rapid Rabi oscillations this well known 

n 

inhomogeneous broadening is suppressed 19]. The slow decay discussed here is due to a 
temporal inhomogeneity, rather than a spatial inhomogeneity in the mean-field energy. The 
origin of this decay is in the spread of dressed states spanning the initial state m z = —j, 
each having a slightly different Rabi frequency. The decay can be explained in the Bloch 
sphere picture by a spread in the initial conditions due to quantum uncertainty. The initial 
state m z = —j is classically described as a vector pointing down on a Bloch sphere of radius 
j. However this is only the expectation value of the angular momentum, and there is a 
fluctuation (independent of j) in the initial value of j x and j y . This uncertainty in the 
initial conditions exists in the linear Rabi oscillation problem too, but all vectors have the 
same oscillation frequency, thus all possible vectors return to their initial value together 
and there is no dephasing. Non-linearity causes the different trajectories to have a slightly 
different Rabi frequency, thus causing dephasing. This decay is therefore not predicted by 
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FIG. 8: The transformation matrices of H2W mt for different values of (j,/Kl = 1,1.5,2 and 4 (a-d 
respectively). For /x < Ml the perturbative classical equation m x = ±\/ j 2 — m 2 + (/x/ Ml)m 2 / Aj 
describes the region of high amplitude in the transformation matrix quite well (dashed line). This 
approximation breaks down as the instability comes in. At /x = 2MI, the turning point reaches 
the initial dressed state vector, when starting from the n = Fock state, and the dynamics depart 
from the perturbative expectation (see Fig. E3) 



the homogeneous Gross-Pitaevskii equation, which corresponds to a single trajectory of the 
expectation value of j on the Bloch sphere. 

The perturbative treatment discussed above is correct for /j < Ml. It is known that for 
/i = Ml a dynamical instability occurs in the solution of the two mode Gross-Pitaevskii 
equation In our model the signature of the instability is manifested by the deviation of 
the region of high amplitude in the transformation matrix from the oval shape, as seen in Fig. 
IH1 For /i = Ml the perturbative classical fit to the shape of the transformation matrix is still 
good, and the curvature of both vanish at m x = j (Fig. |Hk)- For /i > hQ the transformation 
matrix becomes concave (Fig. |Hb)- The right hand side of the transformation matrix varies 
greatly from the corresponding perturbative result, while the left hand side is seen to follow 
the perturbative shape even for /i > HQ. 

We find numerically that for m x larger than a critical m c xl the energy is doubly degenerate, 
with one corresponding eigenvector located only in the top half of the matrix (n < j or 
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m z < 0), and the other located in the bottom of the matrix (n > j or m z > 0). The value of 
m x is found to decrease towards the turning point in the transformation matrix (arrows in 
Fig. |SJ>d) as the matrix size N is increased. Thus we conclude that in the thermodynamic 
limit m c x coincides with the turning point in the shape of the matrix. 

To the left of the instability, eigenvectors have nonvanishing amplitudes both in the 
upper and the lower half of the transformation matrix. An initial state that spans a number 
of such eigenvectors is able to perform Rabi oscillations. The accumulated phases cause 
alternating constructive and destructive interferences for a large range of n. To the right 
of the instability this is not the case. In the thermodynamic limit, each eigenvector is 
located completely either in the top half or in the bottom half of the matrix. Therefore all 
eigenvectors spanned by a state in the upper half will always interfere in the upper half and 
n > j will never be reached. Thus in this region inversion is impossible, an effect known as 
macroscopic self-trapping (31 1. 

The Fock state n = is approximated by a Gaussian superposition of dressed states of 
width yfj, centered at (fi/hfl)j/4. For \i = 2M1, the point of instability is found to reach 
/ij/4. Thus as /i approaches 2Q, the point of dynamic instability reaches the trajectory of 
the initial Fock state n = 0. Fig. El shows the time dynamics of the initial condition n = 
for the same ratios of fi/HQ as in Fig. |HJ For // < 2hQ, we see Rabi oscillations. Even 
though there is a turning point (in Fig. (Hb), as long as it is not on the eigenstates occupied 
by n = 0, oscillations are regular. For \i « 2hQ (Fig. the oscillations decay rapidly, 
and there is a noisy fluctuation with no clear frequency in the population of n. Finally, for 
/i > 2hQ we see the macroscopic self-trapping (Fig. 

In the next section we move on to describe three-wave mixing with the tools developed in 
this section. The three-wave Hamiltonian derived in section lll Al describes interactions also in 
the low k limit, but are limited to weak excitations which do not deplete the macroscopically 
occupied mode. 

IV. THREE- WAVE MIXING 

In this section we discuss three-wave mixing of Bogoliubov quasiparticles, as governed 
by the Hamiltonian in Eq. (|13|) . In the limit of large momentum, where the quasiparti- 
cle operators 6k coincide with the atomic operators a^, the three-wave mixing described 
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FIG. 9: The expectation value for population of state ke as a function of time for the same 
parameters as in fig El For \i < 2KI the dynamics are described quite well by the perturbative 
treatment (a,b). As [i approaches 2MI, the turning point in Fig. IBfc) reaches the initial dressed 
states and the dynamics become unstable (c). For fi > 2HO, (d), we find weak oscillations which 
cannot reach inversion (macroscopic self-trapping). 

here corresponds to four- wave mixing of atomic matter waves Q|. The fourth wave is the 
condensate which is the vacuum of the Bogoliubov quasiparticles. In the limit of low mo- 
mentum, the four- wave picture is incorrect, as each quasiparticle mode &k is a superposition 
of atomic operators with momenta ftk and —Kk. This Hamiltonian does not conserve the 
total number of excitations + n q + nk_ q , indicating these excitations are not particles, 
but rather quasiparticles whose number need not be conserved. There are two independent 
conserved quantities M = n q — nk_ q , and N = n-^ + rik_ q , which lead to the possibility of 
diagonalizing a subspace of Eq. (fT^j) with only one quantum number n = rik- q representing 
the Fock state \(N — n)k, (M + n) q , nk- q ), where n varies between and N. The fact that 
the difference in the occupation of modes q and k — q is constant can be used as a source 
of number-squeezed states. In the low momentum regime suggested here, dephasing is slow 
and pair correlated atomic beams should be easier to generate than in the high momentum 
regime studied experimentally in Q]. 
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A. Resonant three-wave mixing 

We consider first the resonant case A = of Eq. pip, where the sum of energies of the 
quasiparticle modes q and k — q is equal to the energy of mode k, i.e. e& = e q + £|k-q|- 
The modes q and k — q are thus on the energy-momentum conserving surface for Beliaev 
damping from mode k as shown in Fig. EK- H'sw couples between pairs of modes in this 
subspace which have a difference of 1 in n, and can be represented by the (N + 1) x (N + 1) 
tri-diagonal matrix 



' y/{N-0)M ••• » 



K 

Tjres 



y/{N - 0)M y/(N- 1)2(M + 1) 

V( N ~ 1)2(M + 1) 

1 ••• ■•• y/N(M + N) 



\ ^N(M + N) / 

(30) 

with a coupling term K = gy/NoA^/V = fiA ki<i /y/I%. When is diagonalized, we get 
a new set of N + l eigenstates, dressed by interactions, \m x ), where m x varies between —N/2 
and N/2. 

The square of the transformation matrix elements are shown in Fig. E3for N = 100, M = 
50 along with the eigen-energies. As can be seen by comparison with Fig. El the spectrum 
and the transformation matrix seem similar to those of H^. In fact, as the seed M is taken 



to be larger than N, the factor yM + n becomes more uniform, and can be taken out of 
the matrix (jHOj) as a multiplicative factor. Thus in the limit M ^> N, both transformation 
matrix and eigenvalues are well approximated by the two-wave mixing solution described in 
section EH with Ml w Ky/M. 

Even when M is smaller than N (but still large enough to seed mixing), the factors 
originating from the k — q amplitudes in adjacent matrix elements of H 3 \y, are almost the 
same, y/M + n w y/M + n + 1. We use this to approximate Eq. (J30j) locally as a two-wave 



mixing Hamiltonian multiplied by a local yM + n. 

The transformation matrix in Fig. flW a) is indeed very well approximated by stretching 
the n'th line of the rotation matrix (Fig. Of a)) by a factor \J M + nj\J M + n. n is the n 
for which the derivative of y/N — n\/ M + risjn with respect to n vanishes. For M>iV the 
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FIG. 10: (a) Absolute value squared of transformation matrix of H3W res (Eq. for N = 

100, M = 50 and A = 0. The thin lines are at the rows n = 7 and N — 7. The dashed line 
is the line m 2 + m 2 z x (M + n) (M + n) = j 2 .(b) The exact spectrum (solid line) of linearly 
approximated by y/M + n m x (dashed line). Energy is in units of K. The second order perturbation 
approximation (Eq. I33|) is indistinguishable from the solid line. 
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FIG. 11: The n = 7 (filled circles) and the n = N — 7 (empty circles) rows of the transformation 
matrix presented in Fig. HOf a). divided by \JM + n. The x axis is the energy as calculated by 
Eq. j32J), and divided by yM + n. The amplitudes of the different Fock states are normalized by 
[(M + n)/(M + T1)] 1 / 4 . The line is the same function (Eq. I25|) as in the two-wave mixing in Fig. ^ 

approximation of H^w = ^2W becomes exact. 

As can be seen in Fig. ^2 f° r small n, both rows n and N — n agree remarkably well 
with the function ()25|) when plotted versus the spectrum of H^, as in the case of two- wave 
mixing, even for M < N. 

The spectrum is nearly linear centered around with a slope dE = K\J M + n. This leads 
to oscillations, with oscillation period of t = h/dE in analogy to the two-wave mixing case. 
These oscillations decay due to the fact that the spectrum is not strictly linear. In order 
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to model analytically the nonlinearity in the spectrum, we replace both operators 6 q and fej. 



by the operator y M + &k_ q &k-q where the square root of an operator is determined by a 
Taylor expansion. Since in the Hamiltonian (fTTJ) 6 q commutes with fej^k-q, we approximate 
Eq. in a symmetric form: 

We now define angular momentum operators in complete analogy to Eqs. (jl9H23j) with the 
substitution a kB — > 6 k _ q and a — > 6k- In the angular momentum representation Eq. (j3*Tll 
becomes 

Hw ~ f ( V^+7^X i + + J" VWT+Xj • (32) 



To the lowest order in J z / (M + j) we have ~ which is simply iJffp m Eq. 

where = Ky/M + J. We use perturbation theory to second order in J Z /(M + j), to find 
the corrections to the energy due to the non-linearity. 



3(3 + 1 ) \ , e m : 



fRj + 5- 



(33) 



16(M + j) 2 y 16(M + j) 2 _ 

The perturbative spectrum obtained in Eq. ()33)1 is indistinguishable from the solid line in 
Fig. EHb). 

As in the nonlinear case studied in section llll CI the nonlinearity in the spectrum leads to 
dephasing in the dynamics. Fig. IT2lcompares between the numerical (a) and the perturbative 
(b) results of the population of mode k — q vs. time, with the initial condition rik- q = 0. 
It is seen that the perturbative solution captures the non-linear dynamics quite well even 
for M < N. The rapid oscillation at frequency dE/h are unresolved in Fig. EE but 
the envelope is seen to be in good agreement. As expected from wave mixing, there are 
oscillations between the population of the quasiparticle modes k and k — q. The seed in 
mode q oscillates as well and has the occupation M + n. The initial Fock state is a gaussian 



superposition of dressed states of width yjM] (M + n), thus the decay time is approximated 
analytically by tdecay = h/E ^ j M ^ M+n ^ ~ 32 ^^ in agreement with Fig. H21. There is a revival 
after time t reviva i oc N, due to the finite number of quasiparticles. This revival time however 
is much longer than the decoherence time for realistic values of N. Note that as opposed to 
the interacting Rabi oscillating atoms of section llll C[ in which oscillation period and decay 
time are independent of N, both times obtained from Eq. (}3*3*|) scale as yfj, when the ratio 
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FIG. 12: The expectation value of the occupation of mode k — q as a function of time obtained by 
numerical diagonalization (a) and perturbative calculation (b) of H^yy, for the same parameters 
as in Fig. I1U1 and starting with the Fock state n = 0. Time is in units of h/K. 

M/2j is kept constant. Solving for a large system of size j, would require the diagonalization 
of matrices of size 2j + 1. An approximate solution therefore can be attained by solving for 
moderate f and scaling the result by y/j/j'- 

B. Detuning 

In the previous subsection, the seed, mode q, was on the energy-momentum conserving 
surface of mode k. Exact energy conservation is not a necessary condition for wave mixing 
as also a populated mode with a small detuning A from this surface, as shown in Fig. |2t> 
will lead to wave mixing. Figure ITBI shows the dressed state spectrum obtained by numerical 
diagonalization of (Eq. ([13))) at different detunings for N = 10 and M — 5. 

When a detuning term is added to H 3 w, Eq. (|3^jl is modified to 

H 3W « y (^\/M + j + J z J + + J~ \] M + j + j)j + hA (j + J,) . (34) 

A perturbative approximation to the energy spectrum is obtained in a similar fashion to 
Eq. (HUH), with the Hamiltonian ()34|) and is indistinguishable from the solid lines in Fig. 
IT3l At detunings much larger than K/h, the Hamiltonian (J34*|) is nearly diagonalized in 
the Fock basis n or m z = n — j. The energy spacing, and therefore the energy span of 
the spectrum, increases linearly with the detuning. At large negative detunings, the lowest 
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FIG. 13: The spectrum of H^w for N = 10 and M = 5 vs. detuning. Detuning is in units of Si, 
and energy in units of hVt = K\JM + j. The perturbative result is indistinguishable from the exact 
calculation (solid line). Dashed line is E = NhA, the lowest (largest) energy for large negative 
(positive) detuning. Dotted line is E = 0. For large negative detuning, the lowest dressed state 
corresponds to the initial state n = N, while for large positive detuning, the lowest dressed state 
corresponds to n = 0. By adiabatically sweeping the detuning, one can transfer all population 
from modes k and k — q to modes q and k — q. The inset is the difference between adjacent eigen- 
energies. Non-linearity initially increases with |A|, reaches a maximum which depends weakly on 
M/N and is always at |A| ~ Cl, and then decreases for large |A|. 

energy dressed state corresponds to the Fock state where all of the excitations are in modes 
k and q, whereas at large positive detunings the lowest energy dressed state corresponds 
to the Fock state where all of the excitations are in modes q and k — q. In analogy to 
rapid adiabatic passage between two atomic levels, using a sweep in the laser frequency, an 
adiabatic transfer of the excitations population between the excitation modes of a BEC is 
possible using a sweep in the detuning. Given an initial state where all of the excitations in 
modes k and q at a large negative detuning, an adiabatic sweep of the detuning to a large 
positive value, would adiabatically transfer all of the excitations into modes q and k — q. A 
possible way to sweep the detuning during the experiment would be the use of a magnetic 
Feshbach resonance to tune g j^ . 

In analogy with Eq. ([16)1. we find that to zero order in J Z /(M + j), H^w m Eq. (I34|) is 
given by = HVn 2 + A 2 J h where the axis I points in an angle 9 = 7r + arctan(f2/A) from 
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the z axis in the x — z plane. For nonzero detuning the first order perturbative correction 
to the dressed states' energies is proportional to J Z /(M + j). Therefore as the detuning is 
increased from zero, the spread in the differences between adjacent dressed states energies 
(non-linearity in the spectrum) is expected to grow. On the contrary, for large detuning 
Eq. is dominated by the last term and is proportional to |A|, so the spectrum 
approaches the linear spectrum of two-wave mixing with large detuning. This is simply the 
zeeman shift in the angular momentum representation. We therefore expect there will be an 
intermediate value of detuning for which the spread in energy differences will be maximal. 
The inset of Fig. H31 shows the energy difference between adjacent dressed states for the 
same parameters as in the main figure. The non-linearity of the spectrum indeed increases 
initially with the detuning, due to the quadratic nonlinearity which is added in the presence 
of detuning, and decreases towards A for large |A|. 

The spread in energy differences is manifested as dephasing in the dynamics of the system. 
Figure ITU shows ^k_ q ^k-q^) as a function of time, for three different detuning values, A = 
(dotted line), A = Q (dashed line) and A = 3Q (solid line), beginning from the initial state 
|rik = iV, n q = M, rik_ q = 0) with N = 100 and M = 50. The oscillation frequency is found 
to increase with detuning to a value of fi e // = + A 2 . The oscillations amplitude is seen 
to decrease with larger detuning as 1/A 2 . The decay time of oscillations initially increases, 
and then decreases, in agrement with the spread in the energy differences plotted in the inset 
of Fig. for which the maximal spread is for A/Q « 1.5. We use perturbation theory, in 
analogy to Eq. (jS2J), to calculate the spectrum and t decay in the presence of detuning. These 
perturbative results are found to fit the exact solution very well. 

V. BELIAEV DAMPING 

Thus far we have discussed the time evolution of the system within the dressed state 
subspace, which is defined by the initial population of the k and q modes, and have neglected 
damping into empty modes. In analogy to the treatment of spontaneous photon scattering 
as a transfer between dressed state subspaces of an atom interacting with a laser field [l6], we 
now consider scattering into empty modes as transfer between dressed state subspaces. The 
damping of excitations from the mode k is no longer elastic, but rather carries the energy 
difference between the dressed states among which it occurred. The extra energy is taken 
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FIG. 14: The expectation value of the occupation of mode k as a function of time, for three 
different detuning values, A = 0, A = ft and A = 3f2 (dotted, dashed and solid lines respectively) 
for N = 100 and M = 50. The oscillation decay time for A = Cl is much shorter than for A = 
and A = 3f2, in agreement with the energy differences in the dressed state spectra for the these 
detunings (see inset of Fig. I13|) . Time is in units of 1/0. 

from the interaction between the wave packets in the k and q modes and was initially put 
into the system by the lasers which excited the condensate into the (N, M) subspace. We 
consider only damping from the momentum mode k since the damping rate from momentum 
modes q and k — q is typically much smaller . 

The damping is an irreversible process due to the quasi-continuum of N c modes q' (other 
than q and k — q) on the energy-momentum conserving surface. The transition amplitudes 
between all possible pairs of dressed states under 



reveal the spectral structure of the damping process. 

The decay from mode k decreases N by 1 and does not change the value of M, Hence in 
the angular momentum representation, a dressed state \j,m x ) can decay into \j — 1/2, m' x ). 
The matrix elements of the Hamiltonian Humping between numerically obtained eigenstates 
of H 3 w (Eq. H3J) (j — ^ / 2, m' x \ H damping \j,fn x ), are plotted in Fig. fTBTa) as a function of 
the energy difference between the states E mx — E m ' x for iV = M = 10, A = 0. Transi- 
tion amplitudes are large only between dressed states of neighboring energies in the two 




(35) 
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subspaces. The damping process is schematically shown in Fig. fTBT b). Starting from the 
N = M = 10 subspace with 11 dressed states, Beliaev damping of an excitation from the 
k momentum mode into two empty modes, q' and k — q', will transfer the system into 
the N' = 9, M' = 10 subspace with 10 dressed states. Since the energy spectrum in both 
subspaces is nearly linear, and there is one less state in the smaller subspace, the energy 
differences between states coupled by Hdamping are almost the same for all initial dressed 
states. The decay of each dressed state into only two neighboring states in Fig. lT5Tb) results 
in a structure of a doublet in the Beliaev damping spectrum which is plotted in Fig. Hoi for 
experimentally realistic parameters. 

We choose as a model system a condensate of 3 x 10 5 , 87 Rb atoms in the F = 2, rrif = 2 
ground state. The condensate, which is similar to the experimental parameters of |33j], is 
homogeneous and has a density of 3 x 10 14 atoms/cm 3 . The damping rate per excitation of 
each transition is taken as a Lorentzian with a width and height of T/N^ around E mx —E m i. 
iVk is the average occupation of mode k for the dressed state \m x ). Averaging over all of 
the possible transitions between subspaces, weighed by the initial dressed state occupation, 
Fig. EH shows the damping rate from the N = M = 5 x 10 3 to the N = 5 x 10 3 — 1, 
M = 5 x 10 3 subspace vs. energy, for k£ = 3.2 , q = fc/v2 and A = (solid line). The rates 
are calculated numerically according to the Fermi golden rule jl^ : 

Fm x ,m> x = -J- 2^ I ~~ 1/2; m 'x\ Hdamping \j, ^x)\ 2 5 (efc + E mx — € q — C k _ q — E m i^ . (36) 
q' 

A clear doublet structure is evident. This splitting of the spectrum is analogous to the 
Autler-Townes splitting where the spectrum of spontaneous emission is split due to rapid 
Rabi oscillations [34 1. 

Intuitively we expect a doublet in the spectrum to appear as a result of the temporal 
oscillations in the occupation of mode k. The decay from mode k is modulated, and the 
spectrum of a linear system is a Fourier transform of its time correlation function [35]. 
Temporal modulations, therefore, lead to a splitting of the oscillation frequency in the decay 
spectrum. This argument is strictly true only for a linear system such as H 2 w- 

To understand the doublet structure of the decay spectrum, we use the local similarity of 
H 3W to H 2 w, and approximate the three- wave damping rate by the two- wave damping rate. 
The eigenstates of H 2 w, \j, m x) can be viewed as iV = 2j spin 1/2 quanta, n a = j — m x of 
them pointing in the +x direction (right) and np = j+m x of them pointing in the -x direction 
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FIG. 15: (a) The strength of coupling between pairs of dressed states in the two subspaces 
N = 10, M = 10 and N' = 9, M = 10 as a function of the energy difference between the two states. 
In the linear two-wave approximation (empty circles) the coupling is strictly only to the two states 
closest in energy, and the coupling strength is linear in m x as implied in Eq. I|37|) . For three- wave 
mixing (filled circles), there is a negligible coupling to other states, and the coupling is no longer 
linear in m x .(b) Schematic drawing of Beliaev damping as transfer between the subspaces. 

(left). We use the Schwinger boson mapping (sec. IIII Aj) to associate operators a and $ 
with J x . The state \j,m x ) can be expressed by the quantum numbers \n a ,np), just as the 
state \j,m z ) is equivalent to |rik, 7ik_ q ). The operator which annihilates a spin pointing 
in the +z direction is then given by b\ = 1/ y2(a — p) . We use the Fock representation to 
write (^a — (3 \ \n a , rip) = ^Jn^ \n a — 1, rip) — y/np~ \n a , rip — 1). Transforming this to angular 
momentum quantum numbers using n a = N/2 + m x = j + m x and rip = j — m x , we get 



K \j,m x ) w ^y/j + m x j - ^,m x - - ^ j - m x j - ^,m x + 



(37) 



From Eq. (|37j) it is evident that H damping couples the dressed state \j,m x ) only to the two 
dressed states \j — 1/2, m x ± 1/2). 

We have derived Eq. (pT7j) under the assumption that H 3 w is indeed approximated by 
H 2 w- In Fig- IToTa) we see a comparison between the matrix elements of H damping for 
numerically obtained dressed states \j,m x ) (filled circles) and the two-wave approximation 
(empty circles). The spread in the energy difference is due to the non-linearity in the 
spectrum of H 3W , however H damping is found to significantly couple only dressed states with 
neighboring energies as in the two-wave approximation. 

A perturbative result for the decay spectrum where energies are calculated according to 
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FIG. 16: Damping spectrum between the JV = M = 5 x 10 3 subspace and the N = 5 x 10 3 — 1, 
M = 5 x 10 3 subspace, at A;£=3.2, and q = k/\^2. For these parameters Cl = 2ir x 610 Hz. For 
A = (solid line), the perturbative result is indistinguishable from the numerically calculated 
spectrum. For A = Cl the perturbative result (thin dashed line) is still a reasonable approximation 
to the numerical calculation (thick dashed line). The calculation was done with j = 50 and scaled 
to j = 2500. 

Eq. (JHIjj) and rates according to Eqs. ()3b|37j) is indistinguishable from the numerical result 
in Fig. [Till for resonant wave mixing. 

Instead of one shell, which is the energy conserving surface in momentum space for s-wave 
collisions, splitting in the Beliaev Damping spectrum would direct the colliding atoms into 



two separate shells [12[ . Experimentally, the energy doublet can be observed by computerized 
tomography analysis of time of flight absorption images of the 3WM system j^. 

In the presence of detuning, the Beliaev decay rate is calculated in a similar manner, 
again using Eq. (|36|h A damping spectrum for A = Cl is shown as a dashed line in Fig. 
Ibo! The distance between the peaks increases with the detuning to a value of £l e ff ~ 
\/VL 2 + A 2 . As the detuning is increased, one peak grows while the other shrinks. For a 
positive (negative) detuning, the negative (positive) frequency peak becomes larger, whereas 
the magnitude of the positive (negative) energy peak decreases, for a large detuning, as 1/A 2 . 
The perturbative result (thin dashed line) still fits the numerical prediction (thick dashed 
line) well. Deviations of the perturbative result are due to the non-linearity in m x of the 
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FIG. 17: The decay spectrum for different values of detuning. Darker areas correspond to higher 
transition rates. For large positive and negative detuning, the spectra asymptotically approach 
E = and E = hA ( black dashed lines). The gray dotted line is the perturbative result of Eq. 
(|14l) . The thin vertical lines are at the location of the spectra in Fig. 1161 

transition probabilities between three-wave dressed states, while our model assumes these 
probabilities are linear as implied by Eq. (p?T|) . 

In Fig. El we plot the spectrum vs. the detuning, where energy and detuning are in 
units of hQ. Darker areas in the graph correspond to higher transition rates. For a positive 
(negative) detuning, the center of the positive (negative) energy peak is approximately at 
E = hA (diagonal dashed line in Fig. [T7j) . The center frequency of the negative (positive) 
frequency peak approaches E = (horizontal dashed line in Fig. fTTjl like 1/A. For a 
detuning larger than Q, the center of the line is well approximated by the perturbative 
result given in Eq. ()14j). 

Several additional mechanisms will contribute to the broadening of the two peaks beyond 
what is plotted in Fig. The fact that only the first scattering event occurs between 
the N = M = 5 x 10 3 and the iV = 5 x 10 3 — 1, M — 5 x 10 3 subspaces will broaden the 
resonances. Since the energy splitting scales as \/~N, in an experiment where one scatters 
dN atoms from mode k, this will result in a relative broadening of dN/2N. According to 
the same scaling, an initial coherent, rather than Fock, state will cause a relative broadening 
of \/N /2N. In the laboratory, the finite size and the inhomogeneous density profile of the 
trapped condensate, will cause additional broadening, and often dominate the width of the 
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FIG. 18: The width (filled circles) and the peak location (empty circles) of the spectrum plotted 
in Fig. 1171 in the presence of additional broadening taken to be a gaussian of width f2. Although 
the broadening washed out the splitting, the broadening and resonance shift are still apparent. 

Beliaev damping spectrum. In this case the doublet structure of the resonance may be 
unresolved, however the changes in the transition energies and heights, are translated to 
changes in the resonant frequency and the width of the collisional signal as shown in Fig. 
ITS1 There will be a broadening of the decay spectrum around A = 0, whereas the energy of 
the collisional products will fall on a dispersive curve around the resonant value. Figure EH 
indicates that even though splitting and oscillations may be difficult to observe when other 
decoherence mechanisms dominate over the 3WM coupling, line shifts and broadening can 
be still easily observed under such conditions. 

Another way to measure this spectrum is by probing mode k with a weak Bragg pulse. A 
spectrum obtained by setting the momentum of two laser beams, and sweeping the detuning 
between the beams, around the resonance of a quasiparticle excitation of the macroscopically 
occupied mode k, is similar to that plotted in Fig. When the wave mixing is resonant, 
there is a splitting in the spectrum, and when the wave mixing is non-resonant, one peak 
dominates, and in effect the spectrum is shifted as in the optical ac-Stark shift. We calculated 
the Bragg spectrum for the case where the momentum mode q is probed rather than k, and 
found that the spectrum would be a triplet rather than a doublet. This is in analogy to the 
Mollow splitting in atomic physics jl^l . 37 1. 

In many experiments, the Beliaev damping is the main cause of decoherence. From Eq. 
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(|35|) it is clear that the total damping rate decreases for small values of k. This is due to 
a large reduction in the term A^ q/, and the smaller number of available modes into which 
decay is possible q'. This suggests matter wave mixing experiments should be done in the 
low momentum three-wave regime rather than the high momentum four-wave regime. 

In this section we focused on the Beliaev damping of three-wave mixing. The same 
arguments give rise to a splitting in the energy of Beliaev damping from a rapidly oscillating 



BEC discussed in sections III Bl and IIHI We studied the Beliaev damping of such a system 
experimentally, and found a large deviation from the expected s-wave sphere (38]. 

VI. CONCLUSION 

In conclusion, we use the well known Schwinger boson mapping to describe the spec- 
trum and dynamics of wave mixing. Motivated by matter wave mixing of Bose-Einstein 
condensates, we develop a framework for describing two and three-wave mixing. We first 
present the solution of a noninteracting two-wave mixing problem, corresponding to Rabi 
oscillations between two discrete momentum modes, induced by two-photon Bragg coupling. 
Interactions between modes are then added and the exact solution, obtained by diagonal- 
izing the Hamiltonian is compared to a perturbative expansion. When the Rabi oscillation 
frequency Q is smaller then the chemical potential of the condensate the interactions are 
well described as a perturbation over the non-interacting solution, yielding a simple analytic 
expression for the non-linear quantum dephasing time. This dephasing is due to quantum 
fluctuations and cannot be described by mean field. The spectrum of such an oscillating 
3EC, can be calculated using this framework, and is found to agree with the measured one 



22| . For a larger chemical potential, the perturbative solution breaks down, and differs 



from the exact solution. For Ml = fi the GPE is dynamically unstable, and cannot be used 
to predict dynamics. This case of two-mode dynamical instability is different than that 
of Q, which is an instability of a single Bloch state, and is under current experimental 
investigation. 

We then analyze the spectrum and dynamics of three-wave mixing of Bogoliubov quasi- 
particles over a BEC. We study the spectrum and dynamics by treating the three-wave 
mixing locally as a two-wave mixing problem, multiplied by a factor representing the third 
field. By comparison to direct diagonalization of the three-wave mixing Hamiltonian, we 
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find this approximation to be valid even for small seeds M < N. We derive analytic expres- 
sions for the spectrum and decay time of the wave mixing oscillations at different detunings 
from the energy- conservation condition. We consider the coupling to the empty modes, and 
describe the damping from an oscillating mode due to collisions with the BEC. We use the 
local similarity to two-wave mixing to explain the underlying cause for the splitting in the 
spectrum of the Beliaev damping products obtained numerically from the exact three-wave 
diagonalization. 

The Beliaev Damping phenomena studied in this paper, appear also in two-wave mixing. 



A Bose-Einstein condensate undergoing rapid Rabi oscillations, as presented in section IIIB| 
exhibits splitting in the Bragg spectrum of the different momentum modes, and in the Beliaev 
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damping as well. Both these phenomena have been recently observed experimentally 

This dressed state approach can be extended to four-wave mixing of momentum modes 
k 1 ,k 2 k 3 and k 4 , where k x + k 2 = k 3 + k 4 [1]. Here the conserved quantities are N = 
^1 + ^2 + ^3 + ^4, ni — ri2, and n^—n^. If initially two modes are macroscopically occupied 
(say 1 and 3), the four- wave mixing is again similar to two- wave mixing of modes 2 and 4. 
Hence, one observes oscillatory dynamics which lead to a splitting in the decay spectrum. 

Another case of three- wave mixing of bosonic matter waves is binary molecule formation. 
The Hamiltonian governing this process has been studied analytically for the case where 
both atoms are in the condensate mode j^j]. Our model corresponds precisely to molecules 
formed of two separate atomic modes, and can be used to described the evolution of such a 
system. 

This work was supported in part by the Israel Ministry of Science, the Israel Science 
Foundation and by Minerva foundation. 
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